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Abstract 


Landslides can have significant and pervasive impacts to life and property around the world. Several attempts have 
been made to predict the geographic distribution of landslide activity at continental and global scales. These efforts 
shared common traits such as resolution, modeling approach, and explanatory variables. The lessons learned from 
prior research have been applied to build a new global susceptibility map from existing and previously unavailable 
data. Data on slope, faults, geology, forest loss, and road networks were combined using a heuristic fuzzy approach. 
The map was evaluated with a Global Landslide Catalog developed at the National Aeronautics and Space 
Administration, as well as several local landslide inventories. Comparisons to similar susceptibility maps suggest 
that the subjective methods commonly used at this scale are, for the most part, reproducible. However, comparisons 
of landslide susceptibility across spatial scales must take into account the susceptibility of the local subset relative to 
the larger study area. The new global landslide susceptibility map is intended for use in disaster planning, situational 


awareness, and for incorporation into global decision support systems. 
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1. Introduction 


Landslides cause thousands of fatalities annually (Petley 2012; Kirschbaum et al. 2015b; Haque et al. 2016), as well 
as substantial property damage. The true risk may be higher than that observed in recent landslide catalogs due to 
the fact that most casualties are caused by rare catastrophic events (Petley et al. 2005). The first step in 
characterizing the potential impact of landslides (defined in this paper as any mass movements, including shallow 
debris flows, rock falls, and deep-seated rotational slides) is to identify where these events have occurred in the past. 
An ideal landslide inventory would provide both spatial and temporal information on all previous landslides over a 
certain domain. However, most inventories are limited to a short time period that may not fully reflect the 
probability of catastrophic landslides. In addition, many landslides go unreported. Therefore, it is helpful to consider 
not only the historical record of landslide occurrences, but also account for general principles of slope stability when 


predicting the spatial patterns of future landslide events. 


Small-scale (defined in this paper as less than 1:1,000,000 scale) landslide susceptibility maps suffer from four main 
problems: 1) the lack of comprehensive and unbiased landslide inventories; 2) the coarse resolution or absence of 
data inputs; 3) regional differences in the importance of causative factors; and 4) the dearth of expertise on 
landscape processes across a vast region. This work addresses several of these limitations through a heuristic 


approach to represent relative landslide susceptibility at the global scale. 


There have been several projects to represent susceptibility at continental or global scales (Table 1). The Landslide 
Overview Map of the Conterminous United States (Radbruch-Hall et al. 1982) was produced prior to the widespread 
use of digital elevation models (DEM), an input seen in nearly all later research. The European and Indian maps 
(BMTPC and CDMM 2003; Giinther et al. 2014) represent collective efforts in which multiple local datasets were 
assembled into a continental view of landslide susceptibility. In contrast, China and the Caribbean region were 
analyzed as single units (Liu et al. 2013; Kirschbaum et al. 2015a). No variable or method was adopted by all 
authors. The general trend over time is towards “objective” and away from “subjective” methods. Despite covering 
vast areas, all but one of the maps were developed with reference to a relatively small landslide inventory for 
validation. Given this challenge, it is not surprising that most methodologies rely on heuristic methods. Slope and 
geological classification were used most often, while land cover and seismicity were each used in half the studies. 
Classification and ranking of slope data was fairly consistent amongst the studies. This work adopts the standard of 
practice in previous research at continental and global scales and applies it to previously unavailable datasets with a 
flexible fuzzy method. The resulting landslide susceptibility map forms one component of a global decision support 


system that identifies landslide potential in nearly real-time, in concert with satellite-based precipitation estimates. 


32 Table 1: Summary of selected landslide susceptibility maps. In general, previous researchers used heuristic methods to produce 


33 maps with resolutions of approximately one kilometer (30 arcseconds). Some variables, such as slope and geology, were 
34 commonly used, and represent a rough consensus on the importance of these factors. 
Study area USA USA Europe India China Indonesia Caribbean World World 
Radbruch- BMTPC and 
Brabb et al. Giinther et Liu et al. Cepeda Kirschbaum. Nadim et al. Hong et al. 
Reference Hall et al. CDMM 
1999 al. 2014 2013 2010 et al. 2015 2006 2007 
1982 2003 
Landslide events Unknown 24,000 102,000 10 1,200 97 318 3,000 555 os 
s 
Resolution (m) N/A 950 1,000 N/A 1,000 1,000 1,000 1,000 27,800 e 
Analytical 
hierarchy Weighted Weighted Weighted Weighted 
Expert Expert . Neural : Fuzzy . : 
Methodology process & linear linear linear linear 
opinion opinion ae network om overlay Rexed aes 
Frequency combination combination combination combination 
ratio 
Slope x x x x x x x x 8 
Aspect x 1 
Curvature x 1 
Relief x 1 
Elevation x x 2 
Geology xX x x xX x x 6 
Geomorphological 
rpholog) x 1 
classification 
Soil type or texture x x x 3 
Soil moisture x x 2 
Land cover/NDVI x x x x x 5 
Distance to 
stream/Drainage x x 2) 
density 
K6ppen climate 
PP x 1 
classification 
Seismicity x x x x x 5 
Precipitation x x x 3 
Road presence x 1 


35 2. Data 


36 In previous studies (Table 1), slope, geology, land cover, and tectonic features were used most frequently to develop 
37 most of the small-scale landslide susceptibility maps. The same variables were considered in the current work, but 
38 with different data sources (Table 2). Information on roads was also incorporated due to the association between 

39 roads and increased landslide rates (Larsen and Parks 1997; Petley et al. 2007; Kirschbaum et al. 2015b). Global or 


40 nearly global data is available for all of these variables, often without charge. 
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Table 2: The global landslide susceptibility map was created by combining information from four principal sources of 
information for five explanatory variables: slope, distance to fault, geological classification, presence of roads, and forest loss. 


Data Type | Data Set Resolution bo enAalony Extent Source and Details 
Accuracy Variable 
(de Ferranti 2014a) 
Viewfinder Re ae 84 degrees | derived from 3-arc- 
Elevation Panoramas Digital (~90 m) Slope N - 72 second SRTM DEM 
Elevation Data degrees S and several other 
sources. 
Faults and | Geological Map plone A Fa 
Geologic of the World, 3rd | 1:50,000,000 J lonical Global (Bouysse 2009) 
Regions edition cine eae 
classification 
(OpenStreetMap 
; Presence of contributors 2015) Data 
Roads OpenStreetMap Variable Pade Global aggisiis OSM oT 
June 4", 2015. 
Boece Global Forest 30 80 degrees 
Gover Change 2000-— meters/99.6% Forest Loss N - 60 (Hansen et al. 2013) 
2013 ee degrees S 


2.1 Topography 

There are relatively few sources of topographic information with global coverage. One of the best is the Shuttle 
Radar Topography Mission (SRTM). This dataset was initially released at a 3 arcsecond (approximately 30 meters) 
resolution (Rabus et al. 2003), but has been released recently at a resolution of | arcsecond (approximately 30 
meters) and is available from 60° North to 56° South. Unfortunately, the Middle East was not available at this 
resolution at the time of writing. SRTM data contain substantial voids. Several attempts to address this problem have 
been made, including the SRTM 90m Digital Elevation Database v4.1 (Jarvis et al. 2008), Global Land Survey 
Digital Elevation Model (USGS 2008) and HydroSHEDS (Lehner et al. 2008). While many of the SRTM void- 
filling techniques produce reasonably accurate elevations in flat areas, slope and other elevation derivatives can be 
severely affected—especially in mountainous terrain. Each product was evaluated by calculating slope over test 
areas in the Himalayas and the Sahara (where SRTM voids are common). The best global digital elevation model 
(DEM) for the purpose of calculating slope in complex topography was found to be Viewfinder Panoramas (de 
Ferranti 2014a). This is attributed to the use of several sources of topographic information in addition to SRTM, 


which are described below. 


In order to better represent the size and shape of complex topographic features, de Ferranti (2014) reviewed multiple 
series of topographic maps, as well as data from SRTM, the Advanced Spaceborne Thermal Emission and 
Reflection Radiometer (Hirano et al. 2003), the Ice, Cloud, and land Elevation Satellite (Schutz et al. 2005), and the 
RADARSAT Antarctic Mapping Project (Jezek 2002). Landsat imagery was also consulted. Then these data sources 
were combined in a manner designed to draw on the advantages of each (de Ferranti 2014b). Typically, SRTM 


DEM 1-degree tiles with 3-arcsecond resolution formed the basis for the map. Voids in each tile were filled by the 
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most accurate alternative source. The first step in filling voids was to calculate topographic contours from the SRTM 
DEM. Next, the contours were connected across the no-data regions by referencing topographic maps, including 
spot heights. Then the map was searched for artifacts, which were corrected by hand. Finally, the contours were 
converted back to a raster DEM. In some cases, voids were filled directly with data from the ASTER Global DEM 
(GDEM), then checked for artifacts. In Europe, most elevations are based on topographic maps or more precise 
sources, rather than SRTM data. Unfortunately, the tile-based contouring process seems to have introduced errors 
along some tile edges. The specific reason for this behavior is unclear, and the global effect is relatively minor, but it 
should be noted. Nevertheless, this process produces a global DEM with far better representation of SRTM no-data 


regions than other free elevation datasets. 


DEM tiles from Viewfinder Panoramas were converted to slope with R’s raster package (Hijmans 2015), then slopes 
were aggregated to the output resolution by selecting the maximum slope value from the collection of pixels. 
Maximum slope was chosen to represent the most extreme conditions within each pixel and to ensure that the map 
identifies all possibly susceptible areas within the cell. Finally, all tiles were merged into a single map. This slope 
map resembles the Global Slope Dataset for Estimation of Landslide Occurrence Resulting from Earthquakes 
(Verdin et al. 2007), but with the advantage of increased accuracy at high elevations. In addition, the new slope map 
improved the representation of coastal terrain, which should aid decision-making in the British Isles and other 
locations where coastal bluff collapse is a major hazard (Figure 1). However, comparison of the slope datasets 


revealed less than one degree of difference for most locations. 
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Figure 1: The coastline of the Irish Sea and Scottish Highlands as modeled by A) the maximum values from the 30-arcsecond 
Global Slope Dataset for Estimation of Landslide Occurrence Resulting from Earthquakes and B) the maximum slopes 
aggregated from the Viewfinder Panoramas DEM. C) To determine the relationship between these layers, the Global Slope 
Dataset was subtracted from the new slope map. Purple indicates higher slopes in the new dataset. Brown indicates locations 
where the Global Slope Dataset is steeper. Slope values are nearly identical (off-white) in most locations. The primary difference 
is the inclusion of many coastal pixels in the newer map. This probably results from the application of a restrictive land mask to 
elevation data prior to production of the Global Slope Dataset. Since many landslides in the United Kingdom and around the 
world occur at coastal bluffs, inclusion of these pixels may help to balance the concerns of interior and maritime regions. 


2.2 Geology 


Many previous susceptibility mapping efforts include soil and/or rock types as explanatory variables, since each 
material has a unique strength, permeability, and stress history. Unfortunately, most geotechnical properties are not 
available on a global basis. In order to represent this factor, the Geological Map of the World (GMW) (Bouysse 
2009) was simplified into the five categories described by Nadim et al. (2006) and converted to a raster file with a 
resolution of 30 arcseconds (Table 3). The rating was rescaled between 0 and | for consistency with other model 
inputs. The rationale for these ratings was that younger rocks tend to be less consolidated than older rocks, and for 
any given age, sedimentary rocks tend to be weaker than igneous and metamorphic rocks. Nadim et al. also pointed 
out that even though lava rocks may be strong, volcanic deposits are often made of interbedded weak materials. In 
addition, chemical weathering and alteration often have a strong effect on volcanic materials, leaving landslide- 


prone soils and rocks (Frolova et al.; Reid et al. 2001). 


Table 3: Lithological classification 


Material, Age Rating (Nadim et al.) Rescaled Rating 
Water bodies Null 0.1 
Greenland ice cap Unknown 0.1 
Extrusive volcanic rocks, Archean-Paleozoic 1 0.2 
Endogenous rocks, Archean-Paleozoic 1 0.2 
Old sedimentary rocks, Archean-Paleozoic 2 0.4 
Extrusive volcanic rocks, Paleozoic-Mesozoic 2 0.4 
Endogenous rocks, Paleozoic-Mesozoic 2 0.4 
Sedimentary rocks, Paleozoic-Mesozoic 3 0.6 
Extrusive volcanic rocks, Mesozoic 3 0.6 
Endogenous rocks, Mesozoic-Cenozoic 3 0.6 
Sedimentary rocks, Paleozoic-Mesozoic 4 0.8 
Extrusive volcanic rocks, Mesozoic-Cenozoic 4 0.8 
Extrusive volcanic rocks, Cenozoic 5) 1.0 


2.3 Seismicity 


Seismicity increases landslide hazard by destabilizing the soil and debris on slopes, introducing additional fracturing 


that can allow water to penetrate and more rapidly influence the subsurface, and creating steeper or more marginal 
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slopes as a result of seismic shaking and co-seismically triggered landslides (Keefer 1994; Okamoto et al. 2013). In 
addition, tectonically active areas may be prone to increased erosion, due to jointing, graben formation, volcanism, 
stresses (Scheidegger and Ai 1986), and uplift (Larsen and Montgomery 2012). To describe these effects, vector 
representations of major faults were obtained from the GMW. The distance to these faults was calculated to create a 


proxy for tectonic activity. 


2.4 Forest Loss 


Land use is commonly used to explain patterns in landslide susceptibility (Korup and Stolle 2014). However, the 
association between specific land cover classes and the probability of landslides is challenging to characterize 
globally. While ontological difficulties may be avoided by use of a single global dataset, some error is likely 
introduced by grouping disparate biota into a relatively small number of classes. More importantly, there has not 
been clear consensus from the research community as to how to weight these classes. Most studies assign high 
susceptibility to urban areas and low susceptibility to forested areas, which might reflect the impact of 
anthropogenic disturbances on slope stability but could also reflect a bias towards urban areas in landslide 
inventories. The relationship between landslide initiation and land cover classes is more ambiguous. Empirically 
fitted weights would seem to obviate a research review, but biases in landslide inventories can generate incorrect 
associations between specific land cover classes (Steger et al. 2016b), as well as support a false sense of confidence 
in the resulting model (Steger et al. 2016a). Finally, it should be noted that land cover is a constantly changing 
variable (van Westen et al. 2008). The changes caused by fires, urbanization, etc. are likely to have more predictive 
power than the static land cover class itself. For these reasons, land use/land cover was eschewed in favor of forest 


loss. 


Vegetation contributes to slope stability by binding soil particles together and enhancing evaporation (Sidle et al. 
1985; Haigh et al. 1995; Sidle et al. 2006). In a few cases, vegetation may increase hazard, but most slopes are 
strengthened by the presence of vegetation and weakened by its loss. To represent this variable, a Landsat-based 
global map of forest loss from 2000 to 2013 was evaluated (Hansen et al. 2013). The 30-meter forest loss pixels 
were aggregated to a resolution of 30-arcseconds by treating the binary output pixel as “forest loss” if it contained 
any 30-meter forest loss pixel. The resulting map represents forest cover change due to many causes, including 


timber harvesting, fire, and storms. 


2.5 Roads 


Roads may increase the frequency of mass wasting events (Haigh et al. 1989; Larsen and Parks 1997). Particularly 
in developing countries, roads built into and along steep mountain terrain often serve to destabilize the slope (similar 
to ariver cut at a slope’s toe), which can increase the frequency of landslides. After visual comparisons with VMAP 
Level 0 (NIMA (National Imagery and Mapping Agency) 1993) and gROADS (CIESIN and ITOS 2013), the vector 
dataset OpenStreetMap (OSM) (OpenStreetMap contributors 2015) was selected to represent this factor, due its 
more comprehensive and accurate coverage. This roadway network was converted to a raster layer at a resolution of 


30 arcseconds. Larsen and Parks (1997) observed that landslide scars were far more common within 85 meters of 
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roads. While rates remained slightly elevated at greater distances, the effects beyond 100 meters from the road were 
less pronounced. Researchers working at the local scale typically classify distance to road by tens or hundreds of 
meters when mapping landslide susceptibility (Ayalew and Yamagishi 2005; Weirich and Blesius 2007; Dahal et al. 
2008; Regmi et al. 2013; Bhatt et al. 2013; Rubel and Ahmed 2013). With a pixel size of approximately | square 
kilometer, the current susceptibility map cannot model the effect of road construction with the same specificity as 
local studies. Thus, the raster representation of road-related hazards was simplified to the presence or absence of a 


highway in any given pixel. 


2.6 Landslide Inventories 


Landslide inventories from several different events, geographic regions, and methodologies were obtained for 
validation of the global landslide susceptibility map (Table 4). Of these, only the Global Landslide Catalog (GLC; 
Kirschbaum et al. 2015b) covers the entire study area. The GLC was compiled from media reports, online disaster 
databases, and other sources when available, with an emphasis on rainfall-triggered landslides. The database has 
reports from 2007 to the present. In order to reduce the effects of spatial error on validation statistics, 1,194 rainfall- 
triggered landslides with a spatial accuracy of one kilometer or better were selected from a total of 6,790 events in 
the complete GLC. The remaining points could be useful for evaluating products with a coarser resolution, such as 
landslide susceptibility by state, province or country, but were not used for this analysis. Other inventories were 
selected to represent different geographic areas and compilation methodologies (Table 4). Some inventories are 
quite large (Guzzetti et al. 1994; DOGAMI 2015), some cover a long time period (Devoli et al. 2007a; Gerencia de 
Geologia 2012), and some are relatively unbiased but correspond to a single event or observation period (Bucknam 
et al. 2001; ICIMOD 2010). In every case, the local inventory contains more landslides per square kilometer than the 
GLC, which indicates substantial underreporting at the global scale. Reporting biases and uncertainty in this catalog 


have been described in Kirschbaum et al. (2010; 2015). 


Table 4: Landslide inventories used for validation of the landslide susceptibility map. 


Data Set Number of Geographic Source 


points/polygons Extent 


Landslides Triggered by Hurricane 11,555 landslide Eastern (Bucknam et al. 2001) 

Mitch initiation points Guatemala 

Historical landslides in Nicaragua 19,565 points Nicaragua (Devoli et al. 2007b; Devoli et al. 
2007a) 

Landslide Inventory of El Salvador 129 points El Salvador (Gerencia de Geologia 2012) 

Landslide Maps of Utah 2,120 polygons Utah, USA (Elliott and Harty 2010) 


166 


167 


168 
169 
170 
171 
172 
173 
174 
175 
176 


177 
178 
179 
180 
181 
182 
183 


185 


184 


186 
187 
188 


Statewide Landslide Information 
Database for Oregon, release 3.0 
(SLIDO-3.0) 

AVI 


Badakhshan Province Inventory 


Koshi Inventories 


GLC 


3. Methods 


12,095 points 


12,224 points 


609 polygons 


3,407 polygons 


1,194 points 


Oregon, USA 


Italy 


Badakhshan, 
Afghanistan 
Koshi River 
Basin 


Global 


(DOGAMI 2015) 


(Guzzetti et al. 1994) 


(Zhang et al. 2015) 


(ICIMOD 1992; ICIMOD 2010) 


(Kirschbaum et al. 2015b) 


A heuristic fuzzy approach has been taken at the continental (Kirschbaum et al. 2015a), regional (Ahmed et al. 


2014), and local scales (Champati ray et al. 2007), but it has not been previously applied at the global scale. Fuzzy 


landslide models offer some advantages, which include the ability to combine similar datasets in a nested sequence 


prior to the final combination, the ability to use both continuous and discrete inputs, and widespread integration into 


GIS software. A disadvantage is that the output is a “possibility,” which is not strictly comparable to the 


probabilities generated by classical statistics. The heuristic fuzzy approach also enforces transparency, because all of 


the transformation functions are defined in advance. Unlike some machine-learning models, the hypothesis 


represented by the fuzzy overlay model must comport with prior knowledge, not just fit the data. This advantage is 


particularly important for landslide inventories that are known to have significant spatial biases. 


Applying fuzzy logic within a GIS requires two distinct steps (Bonham-Carter 1994). In the first step, a fuzzy 


membership function is assigned for each variable. This function serves to transform the values of the explanatory 


variable to a range between zero and one. The transformation should reflect the relationship between the variable 


and landslide susceptibility. For example, slope was assigned the “large” function available in ArcGIS 10.2 (ESRI 


2013) to represent the fact that susceptibility grows quickly between 10-degree and 30-degree slopes. The second 


step in fuzzy overlay is to combine the fuzzy membership values with a fuzzy operator such as “fuzzy and” or 


“fuzzy or’. In this study, all variables other than slope were combined with the fuzzy gamma operator: 


w=(1- [Jo — fy)” x {Joo 


(1) 


where it is the possibility that a pixel is susceptible to landslides, n is the number of variables to be combined, 1; is 


the possibility that a landslide will occur given the value of the variable i, and y is the parameter that controls 


whether i is closer to the largest or smallest 1;. In order to determine an appropriate value of gamma, it was varied 
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between 0 to | at intervals of 0.1. Each version of the global map was assessed, as described below. Finally, slope 
was introduced through a product function to ensure that no flat terrain would be given high susceptibility values 
(Figure 2). The use of slope gradient as a critical predictor means that the landslide susceptibility map should 
provide more information on mass movements that require a minimum gradient, such as rock falls and debris flows, 


than on low-angle movements such as lateral spreads. 
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Figure 2: The fuzzy overlay model combined data on bedrock, faults, forests, and roads with fuzzy gamma operator, where 
gamma = 0.9. Then slope was introduced with a product function to ensure that no flat ground was identified as highly 
susceptible. Ovals indicate data, while rectangles with the hammer symbol indicate tools from ArcGIS 10.2 Spatial Analyst. 


In order to aid interpretation of the global landslide susceptibility map, the susceptibility values output by the fuzzy 
overlay model were classified into five categories: Very Low, Low, Moderate, High, and Very High. The classes 
were divided at the following fuzzy susceptibility values: 0.11, 0.49, 0.671, 0.75. This classification scheme was 
designed so that each category was twice as large as the next highest, e.g. the Very Low category contains roughly 
twice as many pixels as the Low category. The decreasing category sizes should enable the user to focus efforts 
upon the most susceptible areas. While much of the world is somewhat susceptible, only 3% is very highly 
susceptible. Hazards in this relatively small area may be studied or remediated with greater intensity. Figure 3 shows 


the proportion of GLC locations in each susceptibility category. 
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Figure 3: Distribution of susceptibility for the locations recorded in the GLC (red) and for other areas (blue). The classes were 
divided at the following fuzzy susceptibility values: 0.11, 0.49, 0.67, 0.75. 


Receiver operating characteristic (ROC) curves are commonly used to evaluate the performance of binary 
classifiers, i.e. tests that divide inputs into two outcomes (Zweig and Campbell 1993). Since landslide inventories 
are rarely complete, some locations are likely to contain unreported landslides. This is especially true for the current 
study area, where landslides have been recorded in less than 1% of the map’s pixels. Thus, ROC analysis will give 
only a rough guide to map performance, and other aspects of a landslide susceptibility map should also be 
considered. ROC curves were created for each of the landslide inventories described in Table 4 by calculating the 
number of historical landslides predicted by each possible susceptibility threshold (true positive rate) and the 


number of pixels above each susceptibility threshold (false positive rate). 


A preliminary ROC analysis indicated that low gamma values generated a susceptibility map with a better fit to the 
GLC. However, inspection of the resulting maps showed that the low-gamma maps were dominated by the linear 
inputs, faults and roads. In contrast, the high-gamma maps identified broad regions of hazardous terrain. This 
discrepancy between quantitative and qualitative results can be explained by the fact that many events reported in 
the GLC are associated with road closures, leading to a false level of confidence in low-gamma maps that emphasize 
this feature. Because no single factor (other than slope, which was overlaid separately) is necessary for a landslide to 
occur, gamma was assigned a value of 0.9, which is consistent with the high values published in several previous 
studies (Tangestani 2004; Champati ray et al. 2007; Srivastava et al. 2010; Pradhan 2011; Alvala et al. 2013; Ahmed 
et al. 2014). 
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4. Results 


The landslide susceptibility map is intended to enhance situational awareness with a consistent global picture of 
mass movements. The map covers the Earth’s land surface from 56 south to 72 north latitude (Figure 4). Each 
continent evaluated has susceptible areas, but the major mountain chains (Himalayan Arc, Andes, Alps, and Pacific 


Rim) dominate the map. 


Landslide Susceptibility 
— Very Low 
Low 
| Moderate 
HH High 
BE Very High 


Figure 4: Global susceptibility map developed using a fuzzy overlay model. 


In order to assess the map’s performance, area under the ROC curve (AUC) was calculated for each landslide 
inventory. Both classified (5 susceptibility categories) and unclassified (continuous susceptibility values) maps were 
analyzed to identify any loss of information from classification of the model output into discrete bins (Table 5). 
Uncertainty from spatial error in the inventories was not analyzed, but only GLC points with an estimated accuracy 
better than one kilometer were used. The AUC for the GLC was 0.82, which indicates a relatively successful 
classification of terrain (Hosmer and Lemeshow 2005; Begueria 2006). Local performance of the global 
susceptibility map ranged from very good (Nicaragua) to poor (Badakhshan) (Table 5). Local landslide inventories 
are typically produced for landslide hotspots. In the context of a global map, the entire study area may be highly 
susceptible. As a result, a global classification may place nearly all pixels in one category, giving the appearance of 
randomness in the ROC. AUC values calculated for the classified map were no more than 0.03 lower than for the 
unclassified fuzzy product, which suggests that the classification process preserved most of the available 
information. 


Table 5: The performance of the global susceptibility map was analyzed with eight local landslide inventories. Performance 
appears to depend upon the specific events recorded in each inventory, not upon the size, terrain, or climate of the study area. 


Inventory AUC 
Study Area Unclassified Classified 
Global 0.85 0.82 
Badakhshan, Afghanistan 0.61 0.59 
El Salvador 0.7 0.69 
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Eastern Guatemala 0.7 0.69 


Italy 0.66 0.65 

Koshi Basin, Nepal-India-China 0.84 0.82 
Nicaragua 0.85 0.83 

Oregon, USA 0.75 0.74 

Utah, USA 0.82 0.81 


5. Comparison with previous small-scale maps 


The new global landslide susceptibility map resembles previous publications, both in methods and results. Landslide 
hotspots were identified by Nadim et al. (2006) in many of the same locations that the current study finds highly or 
very highly susceptible. However, the new map identifies a much larger portion of the world’s surface as highly 
susceptible than was shown as highly hazardous in the map of landslide and avalanche hotspots. The difference is 
probably due to the use of a classification system that relies upon “approximate annual frequency” in the prior work. 
The current study identifies some additional large areas as hotspots, including the Appalachian Mountains in the 
eastern United States, eastern Brazil, and Madagascar, which were previously classified as “negligible to very low”. 
This difference is important because many landslides, including fatal ones, have occurred in places like West 
Virginia, Minas Gerais, and Orissa. The new map also has much in common with the previous global landslide 
susceptibility map by (Hong et al. 2007), including large hotspots in the Andes, Himalaya, and eastern Brazil. The 
most notable differences are the relatively low susceptibility ratings assigned to Indonesia, the Philippines, and New 
Zealand by the earlier map. The distribution of categories differs between the maps, with more pixels rated 
moderately susceptible in the map by Hong et al., and more pixels rated very low in the newer map. The significance 
of this is that very few areas can be excluded from future analysis on the basis of the older map, whereas the new 
global map can be used to exclude a majority of the Earth’s land surface from more detailed study. The spatial 
distribution of fatal landslides (Petley 2012) mostly confirms the patterns seen in all three global maps. Highly rated 
areas with few fatal landslides, such as the Southern Andes and the Canadian Rockies, tend to be sparsely populated, 


resulting in fewer reported fatalities. 


The landslide susceptibility map of Central America and the Caribbean region (Kirschbaum et al. 2015a) should be 
very similar to the global map, despite slight differences in methods and data inputs. Since the maps were produced 
at the same resolution (thirty arcseconds) and with the same number of categories (five), quantitative comparisons 
were not difficult. First, the susceptibility classes were assigned a numerical value from 1| (very low) to 5 (very 
high). Then the map of the Caribbean was subtracted from the global map. The global susceptibility ratings were 
equal to or greater than the prior ratings in almost all locations (Figure 5). This is attributed to the fact that Central 
America is a landslide hotspot (Nadim et al. 2006; Petley 2012; Kirschbaum et al. 2015b). While any given location 
within the region may be relatively less susceptible than other Central American locations, it may be relatively more 


susceptible than the Earth as a whole. This tendency was also seen in the susceptibility maps associated with the 
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local inventories described above (Table 4). The handful of locations where the Caribbean map indicates higher 


susceptibility are probably due to the use of an older slope database (Verdin et al. 2007). 


60% 60% 
= Global Landslide Susceptibility Map = Global Landslide Susceptibility Map 
50% = Caribbean Map (Kirschbaum et al. 2015) 50% » Pan-European Map (Gunther et al., 2014) 
40% 40% 
30% 30% 
20% 20% 
” a | on i 
Very Low Medium High Very High Very Low Medium High Very High 
Susceptibility Susceptibility 


Figure 5: This figure shows the relationship between the global landslide susceptibility map and the maps of Europe and the 
Caribbean. The global map often classified terrain as one category higher than the Caribbean map, but the majority of sites were 
identical in both maps. Less than 1% of the study area was classified higher by the Caribbean map. The global map exhibited a 
lower rate of agreement with the European map. In this case, the distribution of higher and lower values was more even. In both 
cases, the number of large differences between global and regional susceptibility classes was quite small. 


The new global map can also be compared to the landslide hazard maps of Indonesia (Cepeda 2010), which have 
four hazard classes. The hazard map for precipitation-triggered landslides shows the highest hazard in southwestern 
Sulawesi and western Sumatra, with the lowest hazard in eastern Sumatra and southern Borneo. The hazard map for 
earthquake-triggered landslides shows the highest hazard in western Sumatra, Morotai, and the mountains of Papua, 
with the lowest hazard in eastern Sumatra and Borneo. Although it has been subdivided into seismic and 
meteorological components, this is roughly the same pattern seen in its predecessor, the global map of landslide 
hotspots (Nadim et al. 2006). The new global map identifies the same locations as highly hazardous, but extends the 
high and very high classes over much of Indonesia, including eastern Sulawesi and many smaller islands (Figure 6). 
Interestingly, both maps portray Java as less hazardous than its neighbors, despite the preponderance of reported 
landslides. This is probably due to population biases in the landslide inventories, but it might indicate the influence 


of anthropogenic terrain modification on landslide rates. 
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Figure 6: Fuzzy landslide susceptibility in Indonesia 


Comparison with the European map (Giinther et al. 2014) revealed large regional differences (Figure 7). In 
particular, the European map shows higher susceptibility in Italy, Ireland, and the United Kingdom, while the global 
map shows higher susceptibility in Iberia and the Carpathians. The global map shows higher overall susceptibility, 
which would be expected if Europe were more prone to landslides than the Earth as a whole. Using the European 
map as a benchmark, the deviation by class is Very Low: 80%, Low: -56%, Moderate: -28%, High: 57%, Very 
High: -17% (Figure 5). It is interesting to note that although the methods for defining susceptibility classes were 
quite different, the overall distribution of European land among classes resembles the global map as a whole. The 
same is not true of landslides (Figure 3); the better validation of the European susceptibility model can be attributed 
to the use of better lithological data, larger and more spatially precise landslide inventories, and a more 
homogeneous study area. Nevertheless, the European map largely confirms the output of the global susceptibility 
model. Less than 0.2% of the map showed a difference of 4 classes (e.g.: a complete inversion from very low to very 
high), less than 2% of the map showed a difference of 3 or more classes, and less than 11% of the map showed a 
difference of 2 or more classes. In other words, there was no difference in 45% of the pixels, and the maps differed 


by a single class in another 45% of the pixels. 
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Figure 7: The new global landslide susceptibility map and the European map (Giinther et al., 2014) were assigned integer values 
from I (very low susceptibility) to 5 (very high). The European map was then subtracted from the global map, so positive values 
indicate areas where the global map has higher susceptibility and negative indicates where the European map is higher. 
Numerous differences between the global and European landslide susceptibility maps can be seen. The European map exceeds 
the global map over large portions of Great Britain and Ireland (brown). The global map shows higher landslide susceptibility in 
most of Portugal and Spain (purple). 
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This comparison suggests that maps produced with different methods, data, and scope may show largely similar 
results. However, maps focused on specific landslide hotspots are not directly comparable to broader overviews 


unless a single, rigorous classification method was applied to both maps. 


6. Discussion 


While comparison to previous small-scale maps revealed strong similarities, this global landslide susceptibility map 
improves upon prior maps in four important ways. First, several new or updated datasets have been released in the 
last decade. In the current context, the most important of these is a DEM made with high-quality SRTM void-filling 
techniques. Second, the use of a conservative method for aggregating 90-meter slope values means that all major 
topographic features were considered by this analysis. Third, the use of fuzzy overlay preserves the full information 
content of continuous variables like slope gradient. Fourth, the simple classification scheme will be familiar to users 
of other susceptibility maps, but the uneven pixel distribution should draw the user’s attention to the most critical 


sites. 


Nevertheless, several features of the new map may limit its use. First, the resolution of the map is approximately one 
kilometer, and terrain varies significantly within many pixels. The choice to aggregate slope by computing the 
maximum value means that some pixels may contain a very small area of steep terrain, while the remainder is not 
susceptible to landslides. Second, the use of biased and incomplete landslide inventories to evaluate the 
susceptibility map makes the results more difficult to interpret. Although this susceptibility model (Figure 2) was not 
fitted empirically, landslide inventories informed the prior research on which it was based. Third, the Geological 
Map of the World is only appropriate for use over very large areas. At local and national scales, more detailed 
information is often available, but varies in quality, format, and cost. Fourth, this map models all mass movements 
with the same treatment. The real world is more complex, and factors which drive rock toppling in Canada are not 
the same as those which can cause debris flows in New Guinea. Fifth, this map does not provide an explicit hazard 
level in the form of an annual probability of slope failure. Therefore, it is very likely that landslides will occur at 
some date in all of the very highly susceptible locations, but the size, frequency, and timing of those events are not 
known. These limitations suggest that the global susceptibility map is best used for a few purposes: situational 
awareness of global landslide hotspots and potential occurrence, the development of global decision support 
systems, and prioritization of future landslide research. It is not appropriate for decisions about infrastructure design, 


building code legislation, or local land-use planning. 


Excerpts from this landslide susceptibility map have already been used during the period leading up to potential 
disasters. In one such instance, the approach of Hurricane Madeline toward the Hawaiian Islands triggered a request 
for information on the potential for landslides. Although the global map is not tailored specifically to this location, it 
was still the most relevant and detailed dataset available to decision makers. This map has also been applied as one 
component of a global landslide nowcast system (Kirschbaum and Stanley 2016). The nowcasts are issued at two 
levels, high-hazard and moderate-hazard, which correspond to different classes of the global susceptibility map. 


After considering susceptibility, a 7-day antecedent rainfall index is compared to historical precipitation levels to 
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identify hazardous locations in nearly real time. While this system focuses upon rain, other landslide triggers, such 


as melting snow or recent seismicity, could be considered in similar models. 


7. Conclusions 


This research assessed landslide susceptibility at a resolution of approximately one kilometer with nearly global 
coverage. The map was evaluated with one global landslide catalog and several local to regional landslide 
inventories. The geographic distribution of landslide susceptibility is very similar to that in previous small-scale 


maps, with the most dangerous terrain located around the Pacific Rim and along the Himalayan Mountains. Other 


hotspots can be found in Europe, Africa and the Americas. While this map benefited from several excellent and free 


datasets, further improvements to thematic data, particularly in soil mapping of mountain regions and landslide 


cataloguing, would improve the results of any future work. The global susceptibility map might be improved by 


incorporation of any future in-situ and satellite-based datasets with improved resolution, accuracy, or completeness. 


The map may be useful for long-term risk assessment and disaster response planning, as well as in the development 


of real-time hazard models. 
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